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Abstract 

Based on functional analysis, we propose an algorithm for finite- norm 

solutions of higher-order linear Fuchsian-type ordinary differential equations 

M 

(ODEs) P(x, £)f(x) = with P(x, &) : = ( J] Pm(x){£) m ) by using only 

m=0 

the four arithmetical operations on integers. 

This algorithm is based on a band-diagonal matrix representation of the dif- 
ferential operator P(x, though it is quite different from the usual Galerkin 
methods. This representation is made for the respective CONSs of the input 
Hilbert space 7i and the output Hilbert space n° of P{x,£). This band- 
diagonal matrix enables the construction of a recursive algorithm for solving 
the ODE. However, a solution of the simultaneous linear equations represented 
by this matrix does not necessarily correspond to the true solution of ODE. 
We show that when this solution is an £ 2 sequence, it corresponds to the true 
solution of ODE. We invent a method based on an integer-type algorithm for 
extracting only I 2 components. Further, the concrete choice of Hilbert spaces 
T~L and is also given for our algorithm when p m is a polynomial or a rational 
function with rational coefficients. We check how our algorithm works based 
on several numerical demonstrations related to special functions, where the 
results show that the accuracy of our method is extremely high. 
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1 Introduction 

Linear ordinary differential equations (ODE) of the type 

pip&m ■= (f>^) (£f ) f{x) = (1) 

are very important tools in many fields (physics, engineering etc.). In many useful 
cases, the functions p m (x) are polynomials or rational functions. As is well known, it 
is difficult in general to solve them analytically for higher-order cases although there 
are relatively general methods for second-order equations with low-degree polynomi- 
als p m (x), for which we employ hypergeometric functions (or special functions) [2] 
and power series expansions about nonsingular points or regular singular points [Tj. 
(Practically, instead of analytical methods, many kinds of numerical methods have 
been proposed and used.) The aim of this paper is to obtain solutions f(x) of a linear 
ODE (JTJ in a Hilbert space H of functions on R when the equation is of higher order 
and/or the function p m (x) is of higher degree. Solutions with finite norm are some- 
times very important in quantum mechanics (e.g. wavefunctions of particles bound 
by potentials) [3J, and for transit or temporary phenomena in signal processing and 
circuit theory that are almost localized in the time coordinate in many applications, 
for example. 

In this paper, we propose an integer- type general algorithm for solving these 
ODEs, by choosing function spaces and their basis systems appropriately. This 
method is based on a pair of Hilbert spaces % and Ti® with distinct inner products, 
where the domain of the differential operator P(x, 4-) is a dense subspace of % 
and its range is a subspace of Under appropriate choice of these spaces and 
their basis systems which will be presented in this paper, a differential equation 
can be expressed by band-diagonal-type simultaneous linear equations, and all the 
'matrix elements' are rational-(complex-)valued. Moreover, under the same choice, 
all the basis functions are rational functions. In addition, from the properties of the 
basis functions used, this method has a somewhat similar feature to power series 
expansions about nonsingular points or regular singular points, in the sense that the 
solution can be expanded as linear combinations of the powers of a rational function 
of x with rational coefficients. From another point of view, this method is closely 
related to the Laurent expansion and hyperfunctions in complex analysis and to the 
Fourier series, under some changes of variable. Therefore, this method has a 'semi- 
analytical' character and it can be discussed from the standpoint of mathematical 
analysis, though it is a kind of numerical method. 

Since it is difficult to apply analytical methods to general higher-order linear 
differential equations, various kinds of numerical methods have been proposed. One 
group is based on the discretization of coordinate or on the differences or on the rela- 
tions between adjacent lattice points (Runge-Kutta methods, for example). Another 
group is based on finite-dimensional subspaces of an infinite-dimensional function 
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space, such as the collocation method, the Ritz-Galerkin method and the Petrov- 
Galerkin method [I] [5], for example. In this group, many kinds of finite element 
methods [5] [6] have been proposed and used widely and efficiently in many fields. 
These methods construct subspaces spanned by finite elements with very localized 
compact supports. In addition, this group contains a subgroup which uses subspaces 
spanned by globally smooth basis functions [4J such as the Hermite functions. 

The method to be proposed in this paper is similar to the latter subgroup in the 
sense that it is based on the finite-dimensional subspaces spanned by global smooth 
basis functions. 

The proposed method is different from Ritz-Galerkin method, in that the function 
space is different from H, and is wider than % in the proposed method. Here, 
remember that the function space TL corresponds to the domain of the differential op- 
erator P(x, and the function space does to the range of P(x, -jj). The choice 
of different function spaces H and may be possible even for Petrov-Galerkin 
method. However, the method proposed here is quite different from the usual 'stan- 
dard truncation methods' or 'projection methods' such as the Ritz-Galerkin and 
Petrov-Galerkin methods, in respect of the following point: 

Although the Ritz-Galerkin and Petrov-Galerkin methods are based on the solu- 
tions of simultaneous linear equations with a square matrix truncated within a finite 
dimension, the method proposed here is based on finite-dimensional truncations of 
the exact solutions of the infinite- dimensional simultaneous linear equations. There 
is a possibility that our numerical solution coincides exactly with the orthogonal pro- 
jection of the true solution to the finite-dimensional subspace, whatever its dimension 
may be, and we have already had some numerical examples where this perfect co- 
incidence occurs really f7J. In order to realize this direction, we solve simultaneous 
linear equations with a non-square-type band-diagonal matrix, in which its column 
is larger than its row. Since the solution is not unique because this non-square-type 
band-diagonal matrix has a non-trivial kernel, we have to extract one solution among 
the above solutions. In order to resolve this problem, we extract one solution among 
them using a novel method, which will be explained in the latter part of this intro- 
duction. Further, this matrix elements do not change when the dimension of the 
subspace increases. Hence the proposed method provides a recursive algorithm with 
no round-off errors up to an arbitrary dimension for the vectors of the space of exact 
solutions of the infinite-dimensional simultaneous linear equations. 

The method to be proposed has five advantages other than the integer-type prop- 
erty mentioned above, as follows; (1) especially when the ODE has no singular point 
or when the ODE belongs to the Fuchsian class even if it has singular points, this 
method can determine the structure itself of the function space of solutions in % 
of the differential equation, directly from the numerical results. (2) Another advan- 
tage is that the convergence of the error to is guaranteed as the dimension of the 
subspace tends to infinity, and an upper bound of the error can be given for the finite- 
dimensional case. (3) Moreover, it does not require any calculation of large matrices 
(inverse matrix, eigenvector, etc.) for solving our simultaneous linear equations. (4) 
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Another strong point is that the basis functions of H are smooth sinusoidal-like 
wavepackets with spindle-shaped envelopes, which are suitable for the expansion of 
various kind of 'natural' functions decaying as x — > ±00. In this sense, the basis 
functions contain both global and local information. (5) Another strong point is 
that this method requires a small amount of calculations for obtaining high-accuracy 
solutions. For example, when the coefficients in the expansion of a true solution by 
the basis functions decay exponentially, the amount of calculations required by this 
method is almost proportional asymptotically to the cube of the number of required 
significant digits. 

In this paper, we will show the validity of the band-diagonal matrix representa- 
tion, i.e., we will show that the square-summable solutions of the simultaneous linear 
equations according to the band-diagonal matrix always correspond to true solutions 
of the ODE except at the singular points of the ODE. Especially, when the ODE has 
no singular point, we will show the one-to-one correspondence between the true so- 
lutions in T-L of the differential equation and the square-summable number sequences 
satisfying the simultaneous linear equations represented by the band-diagonal ma- 
trix. The larger part of its proofs is based mainly on functional analysis. Similar 
one-to-one correspondence can be proved for the cases of the Fuchsian class by a 
modification even if the ODE has singular points. 

However, the presented method has a pitfall based on finite-dimensional trunca- 
tions of the exact solutions of the infinite-dimensional simultaneous linear equations. 
This pitfall is due to the non-uniqueness of non-square-summable solutions of the si- 
multaneous linear equations because the number of linearly independent solutions of 
the simultaneous linear equation is not smaller than the bandwidth whereas the num- 
ber of linearly independent solutions of the differential equations is not greater than 
its order M. That is, there are solutions which do not correspond to true solutions 
of the differential equations, we call these solutions extra solutions. In order to re- 
solve this problem, we propose a method to to remove the extra solutions effectively. 
This method is based on quasi-minimization of the ratio between a norm sensitive 
to divergence and another norm insensitive to divergence. This quasi-minimization 
guarantees the convergence to of the error in the numerical results; the accuracy 
of the numerical results is sufficiently high, even for finite dimensions. 

For minimization of the ratio between two quadratic forms, the usual method 
is based on the eigenspace of the matrix A~^BA~?- with the two corresponding 
inner-product-matrices A and B. However, it is difficult to apply this method to 
the above problem, due to round-off errors, because these inner-product matrices 
are usually very close to a singular matrix with rank 1. However, in this paper, we 
propose an alternative integer-type method for quasi-minimization, which does not 
require as much calculation as the usual method. This method is based on a kind of 
quasi-orthogonalization of integer-valued vectors, which is realized by an idea that 
is conceptually between the Gram-Schmidt process and the Euclidean algorithm. 

An integer-type recursive algorithm similar to the proposed one may be applied 
also to the Petrov-Galerkin method, in order to calculate the head and intermediate 
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rows of solution vectors of the system of simultaneous linear equations described by 
a large-dimensional band-diagonal square matrix. However, in this case, we have to 
calculate new linear combinations of solution vectors satisfying the final constraints 
(linear equations in the bottom rows). For this calculation, we can solve another 
system of simultaneous linear equations, which is described by another square matrix 
whose dimension is half a band width. However, this matrix is usually very close to a 
singular matrix of rank 1. Moreover, the elements in the bottom rows of the solution 
vectors are rational numbers whose numerators and denominators are huge integers. 
Therefore, this type of Petrov-Galerkin method requires a much larger amount of 
calculations than the proposed method based on non-square matrix and integer-type 
quasi-orthogonalization. 



Table 1: Difference from Galerkin methods for norm- finite solutions (without initial 
conditions) 



Method 


Ritz- Galerkin 


Petrov-Galerkin 


proposed 


n° and H 


same 


same / different 


different 


Corresponding matrix 
with truncation 


square 

(band-diagonal / general) 


non-square and 
band-diagonal 


Extra solutions 


No 


can be removed 


Solution vector 


Eigenvalue and eigenvector 
of finite-dimensional matrix 


Exact kernel vector 
(infinite-dimensional) 



The contents of the paper are as follows; Section [2] explains an abstract framework 
for a general algorithm, using a pair of Hilbert spaces H and with distinct inner 
products. Subsection 12.11 states the basic conditions for the pair of Hilbert spaces % 
and r n ) . This subsection gives a sufficient condition for band-diagonal matrix repre- 
sentation for the ODE. In Subsection 12. 2\ using this band-diagonal form, we provide 
the basic structure of our recursive algorithm with removal of non-£ 2 -components. 
Section 3 presents the concrete choices of Hilbert spaces and basis systems used in 
our algorithm, and checks that they satisfy the conditions given in Section Sec- 
tion 4 gives proofs of theorems mentioned in Section [2J In Section 5, we give some 
numerical examples related to special functions and show how effectively our algo- 
rithm works, where we are successful to solve ODEs in a very high accuracy with 
a relatively small amount of calculations (approximately proportional to a power 
of the number of required significant digits, empirically). Moreover, we show how 
numerical results are successful even for the cases with singular points. Section 6 
discusses a related topic and further extensions of our algorithm. 



Table 2: Components of our algorithm 





Differential operator 


Top of [2~I1 


M 


Order of P(x, ^) 


Top of |2.1| 


Pm(x) 


Polynomial: m-th order coefficient func. of P(x, 4-) 


Top of 12.11 


U 


Input Hilbert space 


Top of |2.1| 




Output Hubert space 


m r lo 1 I 

Top ot|2.I| 


A P 


Operator "H — >■ "H defined as action of -|M 


Top of |2.1| 


A P 


Closed extension of 


Top of |2.1| 


Bp 


Operator "H — >■ 7^ defined as action of 4-) 


Top of |2.1| 


Bp 


Closed extension of Bp 


Top of |2.1| 




Basis of TL 


CI 


e° 


Basis of H 


CI 


fn 


f n := (f,e n ) n ioTfeU 


Theorem |2. 1 1 


f 


Vector representation of {/ n }^o 


after (0} 


b n 

m 


Matrix element for Bp: b 1 ^ := (Bpe n , e^)^o 


C2 


to 


Bandwidth parameter: b^ = for \m — n\ > £q 


C2 


Jo 


Integer s.t. 7^ for any integer m > jo 


C5 


N 


Dimension of subspace where recursion is executed 


C7 and Algorithm 


K 


Dimension of subspace of final approximate solutions 


C7 and Algorithm 


D 


Dimension of solution space of b^fn = 


after © 


Po 


Integer p ■= j + £ - 1 


after © 


' h 2 ,K 


'Truncated norm' for number sequences 





2 Abstract structure of our algorithm 

2.1 General framework for general linear ordinary differen- 
tial equation 

In this paper, we consider a general linear ordinary differential equation given by the 
differential operator P(x, -|M defined on the space of M-times differentiable functions 
C M (R): 

P(x,£)f = 0, (2) 

where M is the order of the differential operator P(x, 4-). The main purpose of this 
paper is to analyze the structure of the solution space of ([2]) on a given Hilbert space 
H that densely contains C M (1R). For this purpose, we define the operator as the 
action of the differential operator P(x, 4-) with domain 

D(A P ) := {/ G C M (R) C H\P{x, ±)f E U}. 
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Then, the linear operator Ap is given as the closed extension of Ap with respect to 
the graph norm [10]. That is, we treat the structure of the solution space of the 
differential equation: 

A P f = 0. 

The main goal is to construct an integer-type numerical algorithm for finding non- 
zero solutions of the differential equation given by the differential operator Ap when 
the original space H is contained in a larger Hilbert space as a set. Here, the 
larger Hilbert space also densely contains C M (R). For this purpose, we construct 
a band-diagonal matrix representation of the differential operator Ap under certain 
conditions. In order to obtain a band- diagonal matrix representation, we introduce a 
linear operator Bp from a dense subspace of % to H®, which is defined as the closed 
extension of Bp with respect to the graph norm of the operator Bp defined by the 
action of the differential operator P(x, 4-) with the following domain: 

D(B P ) := {/ G C M (R) H H\P{x, ±)f G 

In order to using a band-diagonal structure, we introduce three conditions for 
the quintuplet consisting of the linear differential operator P(x, the Hilbert 
spaces H and H?, and their CONSs {e n }^1 and {e2}^L , which is abbreviated 
to (P(x, 4-), H, {e n }£i , H®, {e^}^L ). These conditions are shown to hold in sev- 
eral examples for P(x, 4-) later. In what follows, (•, -) H o and (•, •)% denote the inner 
products of and H respectively. 

CI For any n, e n belongs to D(Bp). 

C2 There exists an integer £ such that 6^ := {B P e n , e^) w o = when \n — m\ > £ . 

C3 There exists a linear operator Cp with domain D(Cp) from a dense subspace of 
n° to U such that ei G D(C P ) and (B P f, ei) n o = (f, C P ei) n for / G D(B P ). 

Due to Condition C3, the basis belongs to the domain of the adjoint operator 
Bp. Under these conditions, we obtain the following. 

Proposition 2.1 A function f of the kernel of Ap belongs to the kernel of Bp. 

This proposition is immediate from the fact that the domain of Bp includes the 
domain of Ap. 

Theorem 2.1 Assume that the quintuplet (P(x, £),H, {e n }~ =0 , ft , {e£}~ =0 ) sat- 
isfies Conditions CI - C3. For any function f of the kernel of Bp and any m G Z + , 
the £ 2 -sequence {f n := (/, e n )-^}^L satisfies 

m+l 

E b ^ = °> ( 3 ) 

n=max(0, m—lo) 



s 



i.e., belongs to the linear space 



oo 



V 



{/:= {/Jn=o | E = (m G Z + ) } 



n=0 



{/| E %/n = ° (rneZ + )}. 



(4) 



n=max(0,ra-<o) 



The proof of this theorem will be given in Subsection 14.11 of this section. Due to 
Condition C2, the dimension of V is finite. 

However, square-summable number sequences satisfying do not always corre- 
spond to functions in the domain of Ap (hence in the kernel of Ap). When the linear 
ordinary differential equation (EJ) has singular points, i.e., its solution has singular 
points, we denote the set of the singular points by S. In this case, the obtained 
solutions do not necessarily belong to C M (R), but they belongs to C M (IR\ S) under 
some conditions. So, H has to include the space C M (R \ S). In order to guarantee 
that the obtained solutions are true solutions, we require another condition: 



C4 For any sequence {/ n }^l G V (~)£ 2 (Z + ), the sum \ , fn^n converges to a solution 



Therefore, if the above condition holds, any a square-summable kernel vector of 
the band-diagonal matrix 6^ gives a solution of ODE (j2J) in C M (R \S)D7i. In the 
remainder of this paper, we often use this vector representation instead of a number 
sequence, for simplicity. 

In the non-singular case, the following theorem holds. 

Theorem 2.2 Assume two assumptions: (1) the linear ordinary differential equa- 
tion (TJ|) has no singular points. (2) the quintuplet (P(x, 4-), T-L, {e n }^L ' ^ > ' { e n}ri=o) 
satisfies Conditions C1-C4. Then, the map f {(/, e n )^}^L provides a one-to- 
one correspondence between the £ 2 -solutions of |3|) and the solutions in C M (W) fl H 



The proof is directly derived from Theorem 12.11 and Condition C4 itself as follows. 
When a function / in C M (R) C\% satisfies the differential equation fl2]), it belongs to 
D(Ap) (c D(Ap)) and satisfies Apf = 0. Then, Theorem 12. II and Lemma [2J] guar- 
antee that {{f,e n ) n }^ =0 belongs to Vf]£ 2 (Z + ). The reverse argument is immediate 
from Condition C4. 

By means of Theorem I2.2[ under C1-C4, the linear differential equation is re- 
duced to the simultaneous linear equations (J3J) with a 'band-diagonal structure' of 
bandwidth 2£ + 1. That is, under these conditions, the problem of finding the so- 
lutions in C M (M.) PI H of the differential equation P(x, 4-)f = is equivalent to the 
problem of finding vectors in the space V fl £ 2 (Z + ). 



iV 




/ G C M {R \ S) n H of P(x, £)f = as -»■ oo for H-norm. 



of m>- 
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Even if the linear ordinary differential equation fl2]) has singular points, we have 
a modification of Theorem 12.21 if the differential equation (J2J) is Fuchsian, whose 
definition is given as follows. 



Definition 2.1 An ODE P(x,£)f(x) := p m (x)(£) m f(x) = with analyt 



M 

o m (x)(4-) m f(x) = with analytic 

m=0 

functions p m (x) (m = 0, 1, . . . M) is called Fuchsian if all of its singular points are 
regular singular points [1] [11] [12] . In this case, the differential operator P(x, is 
called Fuchsian. 

Further, in the following, a singular point of ODE P(x,4-)f(x) = is called a 
singular point of the differential operator P(x, 4-). 

As is well known [TTJ [T2], Fuchsian ODEs satisfy the following lemmata: 

M 

Lemma 2.1 Assume that a Fuchsian ODE p m (x)(-^) m f(x) = with analytic 

m=0 

functions p m (x) (m = 0, 1, . . . M) has singular points z\, z 2 , ■ ■ ■ zn- Then, the func- 

N 

tions T\(x — z n ) M ~ m m (m — 0, 1, . . . M — 1) are holomorphic at z\, z 2 , ■ ■ ■ Zm. 

Pm{x) 



\M-m Pm{ x ) 

n=l 

M 



Lemma 2.2 Assume that a Fuchsian ODE p m (x)(£) m /(i) = has holomor- 



m=0 



phic coefficient functions p m (x) (m = 0,1,... M) on K. Then, the set S of its 
singular points is given by p^j (0) . 

As a corollary, we obtain the following: 

M 

Corollary 2.1 Assume that a differential operator Q(x, J^) := q m (x)(-^) m with 

analytic functions q m (x) (m = 0, 1, . . . M) has singular points z\, z 2 , . . . z^ and all of 
zero points of the coefficient function qM^x) are zi, z 2 , ■ ■ ■ zn whose multiplicity are 
not smaller than M. Moreover, assume that the ODE Q(x, -4-)f = is Fuchsian. 
Then, there exist holomorphic functions q m (x) (m = 0, 1, . . . M — 1) onM such that 

Ni 

q m {x) = q m ( x ) Y\( x ~ Zn ^ m an< ^ ^ e se ^ ^ °f ^ s s ^ n 9 u ^ ar points is given by g^(0). 

n=0 

We additionally assume the conditions: 

/oo 
f(x)g(x)v(x)dx. 
■oo 



C2 + There exists a positive function tr in C (R) s.t. (f,g}no = / f(x)g(x)v (x)dx. 



oo 
oo 
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When Condition Cl + holds, % always includes the space C M (R\S) because the set 
S has zero measure. Now, Theorem 12.21 can be replaced by the following theorem: 

Theorem 2.3 Let Q(x, 4-) be the Fuchsian differential operator satisfying the con- 
ditions of Corollary \2.1\ When the quintuplet (Q(x, ^), H, {e n }%L , "H°, {e£}£l ) 
satisfies Conditions C1-C4, Cl + , and C2 + , the map f i-> {(/, e n )^}^ provides 
a one-to-one correspondence between the ^-solutions of (Gj) and the solutions in 
C M (K \ S) n U of (TJj concerning B Q . 

Hence, when the given conditions hold, in order to solve the Fuchsian linear ODE, 
it is sufficient to extract the subspace V fl£ 2 (Z + ) as well as in the non-singular case. 

M 

However, a general Fuchsian differential operator P(x,-^) := ^p m (x)(^) m 

m=0 

does not necessarily satisfy the above condition. In this case, Theorem 12.31 can be 
applied in the following way. Assume that all coefficient functions p m (x) (m = 
0,1,... M) are holomorphic and Pm(x) has zero points Zi,...,Zn G S with the 
multiplicity fii,...,fi N , respectively. Then, Lemma [27TI guarantees the inequality 
f-n — M. Then, the differential operator 

n 

is Fuchsian and satifies the conditions of Corollary 12.11 So, we can apply Theorem 
12.31 to the differential operator Q(x, 4-) in stead of P(x, 4-). 

Since Condition C4 is assumed in Theorem I2.3[ it is sufficient to show the fol- 
lowing theorem, which will be shown in Subsection 14.21 That is, the combination of 
C4 and Theorems 12.11 and 12.41 yields Theorem 12.31 

Theorem 2.4 Let Q(x, 4-) be the Fuchsian differential operator satisfying the condi- 
tions of Corollary \2.1[ Assume that the quintuplet (Q(x, 4-),H, {e n }^L , W*, { e n}n^=o) 
satisfies Conditions C1-C3, Cl + ; and C2 + . Then, the solutions in C M {R\S) RH 
of always give functions of the kernel of Bq. 

In this theorem, the condition for the multiplicity of zero points is crucial because 
it is needed to expand the domain D(Bq) sufficiently. 

Remark 2.1 In general, the inner product (•, of TL does not coincide with the 
restriction on % of the inner product (•, -) h q ofW*. Further, e% does not necessarily 
belong to the domain D(A* P ). In order to characterize Condition C3, we consider 
the special case when = H, e n = e%, and {A P e n ,e m ) n o = (A P e m , e n ) n o . Note 
that the condition = Ti implies that A P = B P . In this case, if the operator A P is 
symmetric, Condition C3 holds. In other words, if the operator A P is not symmetric, 
Condition C3 does not necessarily hold. In such a case, if we define another linear 
operator A' P whose domain is the linear expansion of {e n }, its closed extension A' P 
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is symmetric, and the solution function of A' P f = satisfies the simultaneous linear 
equations (TJ|). That is, if a general non-zero solution function of Apf = does not 
belong to the domain of A' p , this solution does not necessarily satisfy Later, in 
Remark \3.1\ we give an example of the latter case. 

The remaining tasks are divided into two parts: The first part concerns the 
general theory for our algorithm for solving a linear ordinary differential equation 
based on several conditions. This part is called general theory part. The second 
part concerns how to apply the above general theory for several wide classes of linear 
ordinary differential equations. This part is called application part. 

General part 

Task 1 (Subsection 12.2ft Giving a recursive algorithm for band-diagonal-type simul- 
taneous linear equations under Conditions C1-C4. This algorithm requires 
additional condition C5, which will be given in Subsection 12.21 The additional 
explanation for this algorithm is given in Subsection 12.31 

Task 2 (Subsection 12.3ft Giving an integer-type algorithm realizing the above algo- 
rithm by adding Condition C6. 

Task 3 (Subsection H?T]) Proof of Theorem EH] Showing that any C M (M) solution 
corresponds to an element of V. 

Task 4 (Subsection 14.2ft Proof of Theorem 12.41 Showing the one-to-one correspon- 
dence between V n £ 2 {Z + ) and C M (R \S)nU with the Fuchsian differential 
operator Q given in (J5J). 

Task 5 (Subsections 14.3ft Showing the convergence of the algorithm given in Task 
1. 

Task 6 (Subsections 14.4ft Showing that all of £ 2 components in V can be extracted 
by the algorithm given in Task 1 with additional condition. 

Application part 

Task 7 (Subsections 13. 1113^4"]) Constructing "H, , and their CONSs satisfying Con- 
ditions C1-C6, Cl + , and C2 + for a differential operator P(x, 4-) with poly- 
nomial coefficient functions p m (x). 

Task 8 (Subsection 13.5ft Explaining how to apply the above method to a differential 
operator R(x, 4-) with rational coefficient functions r m (x). 

In order to apply our method to a Fuchsian linear ODE P(x, 4-)f = with 
polynomial coefficient functions, it is enough to construct H, T-r, and their CONSs 
satisfying Conditions C1-C6, Cl + , and C2 + with the Fuchsian linear differential 
operator Q(x, 4-) given in (jSJ). Since the Fuchsian linear differential operator Q(x, ^) 
has polynomial coefficient functions, we can apply Task 7 with replacing P(x, 4-) by 
Q(x, £). Due to TheoremES] any solution of P(x, £)f(x) = in C M (R\pl { (0))nU 
can be obtained by this method. 
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2.2 Recursive algorithm for band-diagonal-type simultane- 
ous linear equations 

In the next step, we consider the algorithm for £ 2 -solution of the band-diagonal 
simultaneous linear equations fl3]). In this subsection, we briefly describe the structure 
of our algorithm for this problem and explain how to avoid the usual pitfalls of this 
method. 

From C3, the simultaneous linear equations (EJ have a 'band-diagonal structure' 
with bandwidth 2£q + 1. This type of system of simultaneous linear equations can be 
solved easily. The simultaneous linear equations fl3]) with C2 have at least linearly 
independent algebraic solutions. The linearly independent solutions of the solution 
space V defined in (J3J) can be solved recursively when the following condition holds 
for the quintuplet (P(x, {e n }» Q , H\ {e£}~ ). 

C5 There exists an integer j G Z + such that &™ + ^° 7^ for any integer m > j 
(m G Z+). 

In this case, the dimension D of V is equal to that of 

PO 

n P0 ^ = {{fn} P n= I E h ™f" = (m = 0, 1, Jo - 1) }, (6) 

n=0 

where p := j + £ — 1 and the truncation operator IT m is defined by 

(nj) n = (^ (7) 

In the following, for simplicity, we sometimes identify H m f with the corresponding 
(m + l)-dimensional vector. This band-diagonal matrix b 7 ^ is illustrated by Figure 

m 

So, we define D linearly independent sequences = {Fn }^ = o, ■ ■ ■ , = 
{F^y^LQ by the following procedure: the first p elements of all sequences 
by D linearly independent vectors of V. The remaining elements Fn with n > po + 1 
are calculated by the recursion 

1 n— 1 

F (d) = L_ V b m , F {d) (8) 

n-lo m =n-2£ 

because b™_ eo 7^ there. (The first procedure to find a basis system is easy; it is to 
solve a system of finite-dimensional simultaneous linear equations.) The following 
theorem follows directly from the construction of F^\ . . . , F^ D \ Therefore, we obtain 
the following proposition: 

Proposition 2.2 Under C5, the algebraic solution space U can be spanned by ... , F^ . 
That is, any algebraic solution f of |3j] can be obtained by a linear combination of 
the basis sequences fu,...,fw. 
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£ + l 



^ * * * * 

*vL* vL* vL* vL* vL* 

*Y* *^ *^f» 

vl* vL* vL* vL* vL* vL* vL* 

^^f* *-J-» *^ 

•J-» vL* vl* vl* vL» vl* *A* vL» 

»vL* vL* vL* vL* vL* vL* vL* vL* 
r *y* 



7o 



vL# vL* vl* vL» vL» vt* vL» vl* 











* ▼ 



vl* vL» vl* vL* vl* vL» vL# 

+ 1 = + ^ * # ^ ^ * * * 



vL» vL» vl* vL* vL* 



T 



* : general elements 2 £ ^^^^^^^V 



non-zero elements 



Figure 1: Figure of band-diagonal matrix 6™ 



However, here is an important pitfall. From the existence and uniqueness theo- 
rems, there are M linearly independent true solutions in C A/ (1R) when there is no 
singular point, i.e., S = 0. In this case, therefore, in C M (M) Pi H, the number of 
linearly independent solutions is not greater than M, which is smaller than D in a 
typical example given in Section El Even though there exist singular points, due to 
Theorem 12.11 the solutions in £ 2 (Z + ) of the simultaneous linear equations correspond 
to the true solutions in C M (K\ S) DT-L of the ODE. However, when a solution of the 
simultaneous linear equations does not belong to £ 2 (Z + ), it usually does not corre- 
spond to a true solution in C M (M, \ S) D % of the ODE. Therefore, we have to be 
careful to this differentiate. 

In general, the solution / obtained by the above recursion is a linear combination 
of these three kinds of components, and it is not so easy to extract the component 
corresponding to the true solution in C M (R\ S) C\%. In the following, we propose a 
method to extract the £ 2 -component, i.e., one element of the subset V r fl£ 2 (Z + ). For 
these purposes, we choose a bounded bilinear form £l(f, g) on £ 2 (Z + ) x £ 2 (Z + ) (and 
the corresponding quadratic form Q(f) '■— Q{f, /) on £ 2 (Z + ) ) and the integers K 
and N satisfying 



7 g e(p 



, n(/)> 

N > K > j + 



£ 
1. 



\f n \ 



(9) 
(10) 
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and define the ratio and its minimum: 

*SW : = t<*f£V\{0} (11) 



\i 2 ,K 

a^ N ■= ^min a^(f). (12) 
fev\{0} 

This definition is well defined because Conditions C2, C5 and ( |T0|) guarantee the 
relation 



k2)if := llnjf/Hp > for/GV\{0}. (13) 
Similarly, we define 

* ( Z(f) ■= for/GV\{0} (14) 



K,co 



mm ^(/). (15) 



/6V\{0} 



Hence, our solution space Vr- is given with the following condition 

V K C U K ((^%)- 1 [Q,c t rg] v ]) U {0} (16) 

Now, we introduce our algorithm to approximately obtain the £ 2 (Z + ) components 
of V: 

Algorithm 

Step 1 Calculation of basis vectors of n po V: 

Find a basis system {F^}^ =0 , . . . , {Fn D ^}^ =0 for H po V in (jBJ) by Gaussian elim- 
ination, where D is determined by its result. This is easy because po is small. 

Step 2 Recursive calculation of basis vectors of H n V (po + 1 < n < N): 

Iterate the recursion (jSJ) for n = po + l,po + 2, . . . , N, in order to obtain a basis 
system {F^ =0 , {F^}^ for II^K 

Step 3 Removal of components from H^V corresponding to non-£ 2 -components in 
V: 

Find a linear subspace Vr of 11^- (^(o'^ N )^ 1 [0,ca^ N ]j U {0}. This process can 
be done as follows. 

Step 3.1 Find a non-zero vector {GrP}%=o in IIr- ({cr^ N ) ~ 1 {0, co^^M . 



Step 3.2 Find a non-zero vector {Gffl}n=o i n (i&^ N ) 1 [0,ca 
|Lr„ | n=0 > 

Step 3.n Find a non-zero vector {G$}n=o i* 1 (^(< j ^n)~ 1 1®i c < 7 

X^n J n= o, • • • , i^n J n= o 



(«) ■ 
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n < 



(a) , 



n < 



WhenH^ ((a^r^ca^]) n < {G^}**, . . . , {gT 1 ^ ^ is an empty 

set, we stop this process. Here, ±,K denotes the orthogonal space concerning 
the inner product of £ 2 ({0, . . . ,m}). 

When the purpose is to calculate the truncated elements £ 2 solution Hk(V H 
£ 2 (Z + )), the error is evaluated by the norm concerning inner product (x,y)p K := 
(Hrx, TlKy)t 2 - Denoting the the projection to W concerning this inner product by 
Pw,k, we can evaluate the accuracy of our result Vk by 

\\Pvne 2 (z+),KX - x\\P,K H „\ 

sup ^ . (17) 

SeV K \{0} \\%\\£ 2 ,K 

However, the subspace Vk is not uniquely defined, and is chosen with the condition 
(1161) . The accuracy of our algorithm should be evaluated with the worst 
follows: 

\\Pvkx — x\\p K 
sup sup L --— — 

^cn y ((^%)-Mo,c4%])u{o } *ev K \{o} W\p,k 
\\Pv,kx - x\\p,k 

= SUp ij-q; 

\\P v>K x - x\\p K 

= sup — — . (18) 

It can be shown that this value goes to as follows. 

Theorem 2.5 For fixed K, when N goes to infinity, the convergence 

\\Pv,kx - x\\p K 
sup - — -r— — > (19) 



holds. 



Hence, since Hrx is close to x for x G V fl£ 2 (Z + ), the above theorem guarantees 
that our algorithm gives a subspace of V whose elements are close to elements of 

vn£ 2 (z+). 
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Conversely, the above theorem cannot guarantee that all of elements is approx- 
imated by elements of Vk- For this purpose, we have to guarantee that Pv,kVk is 
equal to Uk(V fl £ 2 (Z + )), which is shown as follows. 

Theorem 2.6 Assume that fl(x) = When we choose sufficiently large num- 

bers N and cq, then for any N > Nq and c> cq, we have 



for any choice of Vk ■ 

Therefore, when we choose fl(x) = \\x\\jp and sufficiently large numbers N and c, 
we obtain a subspace Vk C V close to T1k(V fl £ 2 (Z + )). 

^From the above reason, the above choice of fl suits our purpose. However, when 
we choose the quadratic form fl different from ||x||^2, the convergence is improved 
in several specific example. That is, for sake of practical accuracy of the solutions 
obtained, a condition concerning a kind of sensitivity of the quadratic form fl will 
be required later. The details will be given in our paper [7], and we will use such a 
sensitive quadratic form in the numerical examples in Section [5j 

2.3 Realization by an integer-type algorithm 

When the quintuplet (P(x, ^), {e n }^L , H®, {e^}^L ). satisfies the following con- 
dition, the above algorithm can be realized by the following integer-type algorithm 
with a small modification. 

C6 There exists a complex number 7 6 C such that 7 6^ G Q + Qi, (m, n £ Z). 

The most crucial part is Step 3. To execute Step 3, a simple method is to 
calculate a vector / which minimizes the ratio cr^(/). Usually, with the matrices A 



and B defined by (A) tj := (H K F®, U K F^) P and (B)^ := fl(U N F^,U N F^), this 



minimization can be performed exactly by calculating an eigenvector of the matrix 
A~zBA~2 associated with the minimum eigenvalue. However, this normal method 
is quite difficult to apply, because the matrices A and B are usually very close to a 
singular matrix with rank 1 due to the most diverging components in V, and hence 
this usual method is particularly subject to the 'canceling' due to round-off errors. 
In order to avoid this, many calculations are required, if we try to find the optimal 
vector with high accuracy by the usual methods. 

We can avoid so many calculations and the 'canceling' due to round-off errors 
in this minimization by carrying out the following quasi-minimization. As is proved 
by a geometrical discussion of the convex set in a more general framework in [7j, 
by means of the Schwarz inequality, an arbitrary orthogonal basis system for UnV 
with respect to the inner product (•, -)q contains at least one vector / such that 



p v ,kV k = u k (v n£ 2 (z+)) 



(20) 




Hence, we have only to take the basis vector with the 



minimum ratio cr^' N (f) among this basis system. 
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However, in order to take an orthogonal basis system, we need exact orthogonal- 
ization, which requires a large amount of calculations. To avoid this problem, we pro- 
pose an alternative method based on a kind of integer-type quasi-orthogonalization 
of a .D-dimensional 'lattice', where the angles between the final basis vectors are 
not distant by more than ( (<< 1) from being exactly orthogonal, by which / G 

°ifjv) ccr K",iv]) * s guaranteed for the basis vector with minimum ratio a^ K ^ 

(f) 

This method is somewhat similar to the 'lattice reduction problem' [T3] |14j . which 
is well known as an NP-hard problem if we require exact minimization of the lattice. 
However, our alternative method aims at closeness to orthogonality rather than ex- 
act minimization of the lattice, only with a small amount of calculations, by means 
of a quasi-orthogonalization algorithm which does not increase the integers used for 
the numerators and the common denominator of complex rational numbers except 
for special cases with bad final orthogonality [7J. So, we can perform Step 3.1 with 
few calculations. 

When the dimension Dp of V n £ 2 (Z + ) is strictly greater than 1, we need to per- 
form Step 3.2, . . ., Step 3. Dp. In the first process in these steps, we need to cal- 

culate U K ((agjv)- 1 ^, c a^ }) n < {G^}fU, ■ ■ ■ , {Gn i_1) }* „ > ±tK , which requires 
orthogonalization concerning the inner product ( , )p k- We apply the above quasi- 
orthogonalization algorithm to this orthogonalization. If the second vector {Gn }^=o 
is not orthogonal to the first vector {Gn'}f =0 , linear combinations of {Gn }n=o anc ^ 
{Gn^}n=o ma Y n °t belong to Yi K f(o"i^]v) _1 [0> ccr /!yv]) However, as is shown Sec- 
tion 6 of [7], if {Gn }n=a is sufficiently close to a vector orthogonal to the first vector 
, any linear combination of both vector belongs to TIk ((^km) 1 IP> c<j k^n] ) • 
Repeating this procedure up to the Dp times, we can find a linear subspace Vk as a 
subset of ^(cr^Jv)- 1 ^, ca { ^ N ] J . 

Therefore, we can realize the above algorithm by an integer- type algorithm with 
few calculations. 



2.4 Possibility of the estimation of accuracy 

In numerical methods, it is important whether or not we can determine the accuracy 
of numerical results. For our method, we will give an upper bound of the norm of total 
errors in [7] . This error bound is a function only of the norm of the truncation error 
due to the components outside the subspace tv^' = Span(eo, e±, . . . , 6k), and all the 
other parameters for the bound than this truncation error can be calculated using 
only the numerical results without requiring any knowledge of the true solutions. 
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3 Function spaces and basis systems used in our 
method 

In this section, firstly, we explain what spaces are used for H and as well as what 
basis function systems are used in our algorithm in Subsection 13. II when P(x, 4-) = 
12m=oPm( x )(-^) m can be written as a polynomial in x and 4-. So, in the singular 
case, the set 5" of singular points is given by ^^(0), that is, it is equal to the set of 
zero points of pm- Next, we explain that the presented examples satisfy Conditions 
C1-C5, C1+, and C2+ in Subsections O (for CI, C1+, and C2+), O (for C2, 
C5, and C6). 13.31 (for C3), and 13.41 (for C4). In Subsection 13. 5[ we explain how to 
apply our method to the non-polynomial case. 



Table 3: Notation in Section 4 





Input space (concrete choice) 




<m 


and 


m 




Output space (concrete choice) 




(ED) 


and 


d25D 




Integer parameter for input space 




flzrj 


and 


(El 


k° 


Integer parameter for output space 




m 


and 


(J25]) 


so 


Obligatory minimum of difference ko — 


t.0 


after ([25} 




Wavepacket function used for bases 






(EHD 




^k,n 


'Sorting map': unilateral — > bilateral 











3.1 Construction of function spaces and completely orthog- 
onal systems 

In order to introduce the spaces "H and H®, we state two definitions. 

Definition 3.1 Define the inner product (among measurable functions onM.), parametrized 
by k G Z, as 



(f,9)(k)-= / f(x)g(x) (x 2 + l) k dx 
Definition 3.2 Define the function space 

Lf k) (R): = | /: measurable j \f(x)\ 2 (x 2 + l) k dx < oo] (21) 

J — oo 

( = {/: measurable I ll/H ( fe) < oo} 
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Then, 



Lf k) (R) C L 2 (K) (R) if k > k and L 2 (0) (R) = L 2 (R) . (22) 



J (k)W ^ -^(k) W 11 ^ - ^ cUlu -^(0) 
Moreover, obviously. 



T 2 



For the spaces "H and introduced for the definition of B in Section [5J we will use 



U = L 2 (ko) (R) with (■, •)« = (■, •)(*„), (24) 
«0 = Lj^(R) with (-, = (•, •)<*$), (25) 



where k^ < k — s and s := max m (degp m — m). Then, Conditions Cl + and C2 + 
trivially hold. 

Next, we will introduce the basis function systems {e n } and {e^} for these spaces. 
To do this, we need to define the following functions: 



Definition 3.3 Define the function 

(x + i) k+l \x + i 

Then 



^,n(x):= , [ X —.\ ■ (26) 



1p kj n G 1pk,ri(x) = 1pk,-n-k-l{x) and (^ fc ,m , tpk,n)(k) = 7T <W (27) 

The last orthogonal relation is derived easily from calculation of complex integrals 
by the calculus of residues. When k > 0, as is explained in Section 2 of the paper 
[TJ, the wavepackets defined by (126"]) are 'almost-sinusoidally' oscillating wavepackets 
with spindle-shaped envelopes ^^(x)! = (x 2 + 1)~ 2 , and their approximation 
(for || • || L 2) to sinusoidal wavepackets with Gaussian envelopes holds for large k. 

For these functions, we have the following lemma, which yields the basis system 
of our algorithm: 

Lemma 3.1 {^J^ipk,n \ii <EZ} is an orthonormal basis of L 2 k ^(R). 

The orthonormal property has been shown in the last property of (|2"7|) . Therefore, 
the proof of completeness in L 2 ^ (R) suffices. This is proved in Appendix |A] from 
completeness of the Laguerre polynomials, whose details are omitted here, because 
the Fourier transform of i/jq^ can be expressed in terms of the Laguerre polynomial of 
degree n. The completeness of {ipk,h n€ Z} for k ^ can therefore also be derived 
by (USD and flU]) ). 

Here we point out some properties of ip k , a defined in Definition 13.21 which will 
be important later. 
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Theorem 3.1 Any integer n satisfies 



ipk,ii(x) = --{ijJk-i,h{x) (28) 

£ipk,ii(x) = ni/jk+i,a-i(x) - (n + k + l)ijjk+i,n(x)- (30) 
This theorem can be derived directly from Definition 13.31 

These functions are used for the basis systems of "H and as follows: From 
Lemma [37T] the following {e n }^ =0 and {e^}^ =0 are orthonormal basis systems for Ti 
and in (12^1) and (|25|) . respectively, i.e. Condition CI is satisfied: 

^fco,n fco ,n alld e n = yl^kln k<>n ( 31 ) 

with h Kn : = L-^J + (-1)"+^ L^J , 

where [aj denotes the largest integer not greater than a. 

The indices of functions in { ipk ,n | n G Z} are bilaterally expressed, while the in- 
dices of basis functions in {e n }^ =0 are unilaterally expressed, and they are 'matched' 
to one another by the one-to-one mapping defined by rik,n in (13"Tj) . In order to avoid 
confusion between them, in this paper, the integer indices with double dots " denote 
the bilateral ones in Z, in contrast with the unilateral ones (without double dots) in 
Z + . For ii, the order of the above 'sorting of the basis' for {e n } is 

" fcp+l _ 1 fcp + l I 1 fcp + l _ 3 fcp+l I 3 fcp+l _ 5 fcp + l ,5 » £ pvpti 

2 2' 2 ' 2' 2 2' 2 ~ 2' 2 2' 2 ' 2' '" even 

fc , while it is"-^±i, -^±1 - 1, -!so±L + 1, -*2±I _ 2> -*s±l + 2 , 
- - 3, + 3, ..." for odd fc . For {e£}, similarly with fc£ instead of fc . 

The sorting in fl3T|) may seem to be somewhat complicated and tricky. However, it 
is necessary in order to guarantee Conditions C2 and C5 later. As is mentioned 
later, other conditions hold, we can apply the algorithm given in Subsection 12.21 to 



the quintuplet (P, Lg^R), { J 1 - ^ feo , Jn= o, & (R), { J 1 - $ ^ }» ). 

In addition, when P(x,-^) = Y^=oPm{ x ){^) m * s a F ucns i an differential op- 
erator with polynomial coefficient functions, we need to be careful concerning the 
choice of k$ . This is because we treat the differential operator Q(x,-£-) instead of 
P(x, -4-), which is given in (JSJ). Although the differential operator P(x, -jM requires 

to satisfy k^ < k — max m (degp m — m), the differential operator Q(x,4-) re- 
quires kg to satisfy k$ < k — max m (degp m — m) — ^2 n (M — /i n ), where ji n is the 
degree of nth zero point, whose definition is given in Subsection 12. II With the above 
condition, we can apply the algorithm given in Subsection 12.21 to the quintuplet 



(Q,^ ) (»),{VJK^,.}^.^ ) (R),{VJ^,ft^}^)- 
3.2 Check of Conditions C2, C5, and C6 



A recursive use of the relations in Theorem 13.11 results in the following lemma: 
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Lemma 3.2 Let ko,j,m G Z + and k G Z. When k, < ko + m — j, t/ie function 
xJ (iz) m ^ ; ' c o,n( ;r ) can ^ e expressed as a linear combination ofip K: r(x) 
(r = n — m, h — m + 1, ... , n + m + A; — K ) whose coefficients are polynomials of n 
and ko with degree not greater than m. In particular, in this linear combination, the 
coefficients of the 'outermost' terms with ip K ,n- m and ip K ,h+m+k -K are 

(■ \ ko-K-j+m / \ j m / . \ ko-K-j+m / -I \ j m 

"2 J (2) n("-* +i ) flnd (2) (2) ho™ ib* 

The proof is directly derived from Theorem 13. H where we apply ( 1301) m times, 
next (1291) j times and finally (|28|) ko — n—j+m times. Here note that ko — n—j+m > 
from the condition. In order to guarantee Conditions C2 and C5, we can derive the 
following theorem from Lemma 13.21 Its derivation is based simply on combining the 
results for the linear combinations of Lemma IB.H which is somewhat complicated 
and is given in Appendix iBl 

Theorem 3.2 When the coefficient functions p m (x) (m — 0,1, ... , M) are polyno- 
mials, the function P(x, 4-)e n (x) belongs to . The quantity 6^ defined in C2 
satisfies the following conditions (a) -(c): 

(a) : b n m = if \m — n\ > 2M + k — k } . 

(b) : There exists a polynomial A(x) of degree not greater than M such that 

|&ml — A(n) for any m,n G Z + . 

( c ) K-(2M+k -kS) ^0forr>2M + k + max(-^, 0). 

This theorem shows that the above mentioned quintuplet satisfies C2 with £q = 
2M + ko — kg. Hence, the dimension D of V is greater than the degree M of the 
differential operator P because D > lo > 2M. This theorem also guarantees that this 
quintuplet satisfies C5 with £ = 2M + k — k^ and j = max(/cQ , 0) if pm(±«) 7^ 0, 
because (2M + ko) — (2M + ko — k^) = k^. Hence, the recursive algorithm Theorem 
12.21 can be applied when p M (±i) 7^ and v i6l Pm(%) 7^ 0. 

The accidental cases where Piw(z) = or Pm(— i) = can be easily avoided by a 
change of coordinate x — > x + b for appropriate b G R, because Pm( x ) has only M 
roots. More generally, we can use a change of coordinate x — > ax + b for appropriate 
a > 0, b G R which is useful not only for this but also for rapid convergence, 
by 'matching' of the scale and the position of the localization between the basis 
wavepackets and the true solutions. 

Many band-diagonal elements vanish in the matrix (6^). Especially, the equation 
b r n = holds when n < k$ — 1 and r > ko- This fact can be shown as follows. In the 
above expansion of x J (-^) m ipk ,h(x), the terms with ij) K ^ (f < —1) vanish when < 
n < m— 1, and the terms with ip Kt f (f > —k) vanish when — ko — m < n < —ko — 1. 
These properties are derived from -^ipk ,o(x) = —(ko + 1) iftk +i,o(%) (without the 
term n^+i.fi-i) and £^ ,-fco-i( x ) = ~( k o + 1) ^fco+i.-fco-^) (without the term 
— (h + k + 1)^0+1, a) which are special cases of fl30l . Applying the matching (l3~Tj) 
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to the above vanishing property for k > 0, we obtain the following. When n > ko, 
the terms in e^, (n' < k — 1) vanish in this type of expansion which 
is derived from {e n \ n < k — 1} = {ipk ,n I — ko < n < —1} and {e^ | n < k — 1} = 
{^K,n | — k < n < — 1 }. So, we conclude that b r n = for n < k$ — 1 and r > & - 

The calculations of b™, := (B P e n , e^,) w o need the recursive use of the relations in 
Theorem 13.11 in the bilateral expression. Its program can be realized as an integer- 
type program under C6 in a practical algorithm explained in the paper [TJ, where 
relations ( 128]) . ( 129]) and ( 130]) are modularized. For the recursion (JBJ) in Theorem 12.2] 
we have only to know that &™lf Q+1 , ^_^ - Hence, Condition C6 holds 

with this quintuplet. Here we omit the 'sorted version' in the unilateral expression 
of (128]) . (|29j) and (130]) . because it is too complicated to use in a practical program. 



3.3 Check of Condition C3 

In order to check Condition C3, we define the operator Cp by 

M degp m 

(C p g) (x) := E E (-i)^^ 2 + ir k0 &(x J (x 2 + l) k »9(x)) 

m=0 j=0 
deg p m 

with p m (x) := p m ,j and domain 
i=o 

£>(<5 P ) = {/ e C A/ (R) n L^ } (R) | C P f e q ko) (R)}, 

and describe its closed extension by Cp. 

Theorem 3.3 Under fcp < k — sq, the operator Cp and the operator Bp defined by 
the action of P(x, 4-) in Section^ satisfy 

*/ G D(5 P ) and ^GZ, (B P f, ^) {ko _ so) = ( /, C^itf, ' 

Theorem 13.31 guarantees C3 under the choices ( |24l) . ( 125]) . and ( |3T]) even when pm has 
zero points, because D{Cp) is dense in L? feo _ so ^(R). 

Since the proof of this theorem requires many pages, it is given in [8]. Here, we 
explain briefly the basic idea used in the proof. The equality 

( Bp /, ip k (> ft ) ( feo _ So ) = f /, C P if) k o ~ J can be shown by iterative use of the 'inte- 

rb rb 

gration by parts' / p{x)q\x)dx = \p(x)q(x)\ — p'{x) q(x) dx if we can show 

J a J a 

the disappearance of the contribution of the term [p(x) q(x)] = at each step of the 
iteration in the limit as a — > — oo and 6 — > oo. We can show its disappearance under 
the conditions in Theorem 13.3] even when p(x) and q(x) do not converge as x — > ±oo, 
by means of a 'modified kind of smoothing operator' T which 'blurs' the endpoints a 
and b so that (T n p) (x) and (T n qj (a?) may converge to as x ±oo for an integer 
n. 
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Remark 3.1 The inequality k$ < k$ — sq required in Theorem Iff. 31 is essential for 
Condition C3. Remember that, even though the ODE can be represented formally by 
a band-diagonal matrix, this band-diagonal representation is not always valid if C3 
deos not hold. As an example, we consider the ODE 
i d \ 

-(x 2 + 1) — — h [k + l)x — a J f(x) = {k: integer, a: rational constant) with 



dx 

the choice k$ = k = k. In this case, Sq = 1, and hence the inequality k$ < 
ko — Sq is not satisfied. While the solutions f(x) = C 1 fc+1 (j^ry) (C: const) 

belong to L? k JM), their corresponding vectors do not satisfy the simultaneous linear 

m+£o 

equations &m/n — when a is not integer and C ^ 0. Hence, due to the 

n=max(0,m— £q) 

contraposition of Theorem \3.3l Condition C3 is not satisfied. In other words, in this 
case, Ap is not symmetric, which has been explained in Remark \2.1[ 

3.4 Check of Condition C4 

Next, we will show that C4 is satisfied under the choices (J2"ll) . f[2"oT) . and fl3T|) . When 
the coefficient functions p m are polynomial, the set S of singular points is given by 
the set of zero points of pm, i-e., p^O). 

In a general framework of the theory of elliptic differential equations, the following 
fact is already known. Assume that a function / e L 2 (IR) satisfies the conditions 
(J^) m / G L 2 (M) for m = 1, 2..., M—l and belongs to the kernel of the closed extension 
of an elliptic differential operator J2m=o r m{ x ){j^) m satisfying 3 d > r M {x) 
Then, the function / is smooth, i.e., belongs to C M (M. \Pm{0)) flL 2 (R). This fact is 
shown by a generalization to higher order cases of the discussions in •!)] , for example. 

However, there are many functions that do not satisfy these conditions even for 

true C M -solutions of ODEs. For example, the function fix) = — cos(x 3 + x) 

3x 2 + 1 

is a true solution of the ODE 



which is equivalent with 



( (3x 2 + 1) 2 (-^-) 2 + 6x(3x 2 + 1) (4-) ~ 6(3x 2 - l)(3x 2 + l) 2 - (3x 2 + l A 
\ dx dx J 



f(x) = 0. 



This solution / can be written as the form J2n=o /« ■y i^fe ,ri fco , n with {f n }^ =0 G 
£ 2 (Z+). 

Hence, in order to show C4, we should check that the general solution function 
/ belongs to C AI (R\p~ A . I (Q)) DL 2 (R). The following theorem is important for check 
of Condition C4. 
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Theorem 3.4 When the coefficient functions p m (x) (m = 0, 1, ...M) are polynomials 
of x and the linear space defined under the quintuplet (P, L^(R), *Pk ,h ko n }??=o> 

L^_ (R), {y^V'fe* no }«=o)> ^ en / or an 2/ element f E U fl £ 2 (Z + ) ; i/iere exzsfo 
V? G a M (R\p]^(0)) such that 

oo 

^/ne„(a;o)= (P(x,f)^)(x ) = 



n=0 



for Vxq G 



r(0). 



Theorem 13.41 seems to provide Condition C4 directly. However, this theorem only 
guarantees the point-wise convergence while Condition C4 requires the convergence 
in the norm of H. Hence, we need the following lemma. 

Lemma 3.3 If there exists a function cp G C M (M. \p~m (0)) such that 

N 

lim \ f n e n (x) = (p(x) holds for any x G M. with a sequence {/ n }^Lo e ^ 2 (^ + ); 

N— >oo — ' 



n=0 



then ^hnj / n e„) - <p 



n=0 



Proof of Lemma 13.31 

Since {/nj^Lo belongs to £ 2 (Z + ) and {e n }^1 is a CONS of "H, there exists a 



TV 



function / such that lim (S f n e n ) — f 



n=0 



H 



0. Hence, there exists a subsequence 



{N„}^L such that lim N / n e„(s) = /(#) (a.e.). Therefore, from the trigonometric 



n=0 



inequality, |/(a:) -(p(x)\ < lim (\ f n e n (x)) - ip(x) +\Q2f n e n (x))-f(x) 



n=0 

(a.e.). Therefore, ||/ — (p\\<H = 0, and hence 

N N 



n=0 



Jim (V f n e n ) - ip < lim ( ( V / n e„) - / 



n=0 



0. 



n=0 



Therefore, Lemma 13.31 and Theorem 13.41 guarantee Condition C4 even when has 
zero points. 



3.5 Application to the non-polynomial case 

In this subsection, we explain briefly how the method proposed in this paper can 
be extended to a more general case where the coefficient functions in the differential 
operator are not necessarily polynomials but rational functions of x. 
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We can generalize the facts shown in the preceding subsections of this section, for 

M 

differential operators written in the form R{x, ^) = ^^r m (x) (J^) m with rational 

m=0 

functions r m (x) (m = 0,1, M). Multiplying the least common multiple of the 
denominators of r m (x) (m = 0, 1, M), we obtain a differential operator P(x, 4-) = 
Ylm=oPm( x ) (^) m w ith the polynomial coefficient functions. Then, we can apply the 
algorithm given in Subsection IO to the quintuplet (P, L 2 ^(R), {J\ ipk ,n kotn }n=oi 
L^(R), { y^V'fc* a o )n=o)- Since Condition C4 holds, the numerical result / is 

close to a solution of P(x, £)f{x) = in C M {R \ ^/(O)) n L 2 ko) (R), which is a 

solution of R(x, ^)f(x) = 0. 

However, any solution of P(x, £)f(x) = in C M (R\plf(0))r}L 2 kQ) (R) is not nec- 
essarily obtained by our algorithm in general. When the differential operator R(x, ^) 
is Fuchsian, all of R(x, £)f(x) = in C M (M \p M 1 (0)) n L 2 ko) (R) can be approxi- 
mately obtained by our algorithm. This fact can be shown as follows. Any solution 
of R(x, £)f(x) = in C M (R \ p^(0)) n L 2 ko) (R) is a solution of Q(x, £)f(x) = 0, 
where the Q(x, 4-) is given in (J5]) from P(x, 4-). Due to Theorem 12.41 any solution 
of Q(x, £)f{x) = in C m (R\pm(0)) n LL(1) can be approximately obtained by 
our algorithm. So, we obtain the above fact. 

When the ODE R(x, £)f(x) = has no singular points, we have the following 
stronger characterization for the ODE. This condition is equivalent the non-existence 
of no zero points in the coefficient function pm- In this case, all of the solutions of 
ArJ = in L 2 ko j(R) can be approximately obtained by our algorithm. This fact can 
be shown from the following theorem. 

Theorem 3.5 Assume that pm has no zero points. For any k > 0, there exists an 
integer k$ such that following conditions for f G L 2 kQ ^{R) are equivalent. 

(1) A R f = 0withU = L 2 kQ) (R). 

(2) B P f = with H = L 2 JR) and = L 2 .(R). 

(3) The £ 2 -sequence {f n := (f, e n )^}^L belongs to V defined with with the quintuplet 

( P > L tk 0) W , { ^o, * k0 Jn=0 , L\ kt) (R) , { ^0, n kt Jn=o)- 

(4) P(x, £)f(x) =0andfe C M (R) n L 2 ko) (R). 

(5) R(x,£)f(x) = and f eC M (R)nL 2 ko) (R). 

Proof: There exist an integer k\ and a constant c such that [x 2 + l) fcl ( ^|^ ) 2 < c. 

Then, we choose k$ with satisfying the condition k$ < min{fco — niax m (degp m — 
m), k + ki}. 
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The property k$ < k +ki yields (1)=^(2). The property k^ < k —m&x m (degp m — 
m), Theorem 12. 1[ and Condition C4 for the above quintuplet imply (2)=^(3)=^(4). 
Since (4)=>- (5)=>- (1) is trivial, we obtain the desired argument. ■ 



4 Proofs of Theorems given in Section [2 
4.1 Proof of Theorem 12.1 



Now, we prove Theorem I2.1[ with the following definition, as follows: 
Definition 4.1 Define 

H (n) := span(e , e x , ...e n ) and := span(e£, ef , . . . e£) (32) 

wi/i / fl3]) . and / T57]) . and define the orthogonal projectors P n and P n to W- n > and 
T-i^ n \ respectively, with respect to the inner products (•, and (■, respectively. 

Proof of Theorem \2.1\ By definition I4.1[ lim ||P n / — /||^ = for fE.fi, and 
lim ||P n / — /||^o = for / G Ti®. Hence, P n f and P n / weakly converge to / with 

respect to the respective inner products. 

The condition C3 holds for every f n in any function sequence {/„ G D(B)} 
converging to / G D(B) with respect to "H-norm, and the definition of the graph 
norm guarantees that Bf n converges to Bf with respect to the "H^-norm. From 
these facts, the condition C3 holds even for / G D(B)\D(B) , and hence it follows 
that v n, el G D(B*) with B* = C. Hence, (B(P m f ), e%) n o = (P m f, B*el) H , which 
implies 

lim (B{P m f), e^U, = lim (P m f, B*e*) n = (/, B*e$)„ = (Bf, e£> w o. 

m— >-oo in— >oc 

Therefore, any solution / G D(B) of Bf = satisfies lim (B(P m f), e„)^o = 

(Bf, e^)^o = 0. On the other hand, from C2, it is easily shown that 
v m>n + £ , P n B{P m f) = P n B(P n+i J). Since P n e% = e%, 

\im(B(P m f),e%) H o = lim (B(P m f), P n et) H , 

m— >oo m— >oo 

= lim (P n B(P m f), e§ n * = (P n B(P n+io f), e° n } H o. 

m— >-oo 

These facts lead us to (P n B(P n+ e f), e^) n o = 0, which is equivalent to (j3J) in The- 
orem [2H1 from CI and C2, because b 7 ^ = for |m — n\ > Iq. Thus Theorem 12.11 
holds under C1-C3. 
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4.2 Proof of Theorems 12.41 



In order to prove Theorem 12.41 we prepare the following definition and lemmata: 

Definition 4.2 For a nonnegative integer L and a positive real number e, define the 
function 



W. 



L,z,e 



X 



\x — z\ 



with 



( 



w L (x) :-- 



l -J 
u L (l - u) L du 



u L (l-u 



(ifx<0) 

) L (ifO<x<V 
(ifx > 1). 



This definition results in the following lemma directly. 

Lemma 4.1 The function Wl, z ^(x) has the followimg properties: 
(a): W LjZ , e eC L (R). 
(J): v i6l, 0<W L , z , e (x) <1 
(c): v x eR\{z-e,z + e), W L;Z , e {x) = 1 



(d): 



> dx i 



l Wr z Jx) 



L. 



= for m — 0,1, 

(e): 3 iT L , m > such that v e >0, v xG R, \{£) m W L)Z ^x)\ < K L , m e~ m 
for m = 0, 1, . . . , L. 

Note that Kl,™, does not depend on z. Lemma 14.11 and the property of Fuchsian 
class lead us to the following lemma: 

Lemma 4.2 Assume that a differential operator Q(x, 4-) is Fuchsian, its coefficients 
functions q m (x) (m = 0, 1, . . . M) are holomorphic on IR ; a Hilbert space 7-L satisfies 
C1+. Choose a solution f G C M (R\q£(Q)) n H the ODE Q(x,£)f = (for 
x E R\g^/(0)). Then, for any integer L satisfying L > M, there exist positive real 
numbers e, a and a positive constant K such that 



z n —0 
z n +e 

Zn+0 



[ X Z Yl 



\m/ d ' 
' \dx' 



(W L>Zn;e {x)f(x))\ dx < Ke a 



hold for any real number e G (0, e), any zero point z n of qM{x), and m — 0,1, 



M. 
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Proof of Lemma \4.S\ 

Corollary 12.11 implies that the set S of the singular points is given by g^/(0). 
In the following, arrange the zero points so that Z\ < z 2 < ...z Nl , and, for a 
convenience, let z = — oo and z Nl+1 = oo though they are not zero points. Since 
the ODE Q(x, 4-)f = is of the Fuchsian type, as is well known in the theory of 
the power series expansion about regular singular points [lj, the solution f(x) can 



written as ^ ( ) r (log(x - z n )) s g n>S) ±(x), with a holomorphic function g n . 



s,+ 



s=0 

and g n ^-(x) defined in (z n ,z n+1 ) and (z n _i,z n ), respectively, about a zero point z n 
of qM{x). Here r is the exponent which is a root of the indicial polynomial and s is a 
non-negative integer not greater than M — 1. Hence, (-£;) m f(x) can be written as a 

s m min(m— u,s) 

linear combination c Uit) (x-z) r ~ u_t ' (log^-z))* - ^ — ) m ~ u ~ v g n ^ + (x) 

s=0 u=0 v=0 

with coefficients c UjV in C for x G (£ n , z n +i). Moreover, since g nj §,+(x) is holomorphic 
at x = z n (n — 1, 2, . . . Ni), there exist positive real numbers e and J n ,m,s,+ such that 



dx' 



'9 



n,s,+ 



< J, 



for m = 0, 1, ... , M; n = 1, 2, . . . iVi; s = 0, 1, . . . , s. Similarly, there exist positive 
real numbers eo and J n m s - such that 



X G ( Zy l 



6, Z^ 



\ ; ) 9n,s, 

dx 



ix 



for m = 0, 1, . . . , M; n = 1, 2, . . . N x ; s = 0, 1, . . . , s. 

Let p be the real part of the exponent r used in the above-mentioned expan- 
sion of f(x). Since x(logx) s is infinite-times differentiable for x G (0, oo) and 
lim x(logx) s = for s = 0,1,..., s, when t is a nonnegative integer, from the 

above fact, there are positive numbers K and e such that 
v x G (z„ - e, 2„) U (z^Zn + e), v m G {0, 1, . . .m}, 



. X Zr 



< K\x- zJ t+ P-™ 



holds for n — 1, 2, . . . Ni if t is a nonnegative integer. This fact and the properties (a), 
(d) and (e) of Lemma 14.11 lead us to the statement that there are positive numbers 
K m fh and e such that 



1 x G (z n - e, z n ) U (z n , z n + e), v e G (0, e) 
d m 



< 



rh=0 



rh.\X Z'j 



Vb+p—m^m—m 



holds for m = 0, 1, . . . , M and n — 1, 2, . . . iVi. 
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Since p > — ~ (otherwise, Condition Cl + would result in / H), the inequality 
m > m — p — 2 holds. Hence, it is easily shown that the lemma holds with a = p + | 
(> 0), because / belongs to H and Condition Cl + is satisfied. 



Proof of Theorem \2.J\ 

Let Zi, Z2, ■ ■ ■ , Z]sr x (z\ < Zi < . . . < z^) be the zero points of qM{x). For a 
solution / in C M (R\q M \0)) of the ODE Q(x,£)f = 0, with e > 0, define the 
function 



f e (x) := (jl W M,zrJx)) f{x). 



Then, from the properties (b)-(d) of Lemma 14.11 and Condition Cl + , it is easily 
shown that 



lim J|/ e (x)-/(x)||« = 



(33) 



without any complicated problem caused by the fact z n e q^(0). 

On the other hand, Corollary 12. 1\ Lemma [4.21 and Condition C2 + lead us to 



lim 

e-S>+0 



II 

n=l 



0. 



(34) 



HO 



The convergences (1331 and (j3"4"|) imply that / belongs to the domain of Bq and the 
equality Bq/ = holds. Since S = (0), we obtain the desired argument. ■ 



4.3 Proof of Theorem 1231 

The definition of £ 2 (Z + ) implies that 

(^Lr 1 [0,co|L]c^n£ 2 (Z + ). (35) 

In order to show Theorem 12.51 we denote the set of normalized vectors in V and the 
e neighborhood of a normalized vector x in this set concerning the norm || \\p^ K by 
Ok and U e> g. Thus, from ( 1351) it is sufficient to show that for any e > 0, there exists 
an integer Nq such that 

Ok n (agkrMO, cogv] C 17 := |J ^ 

for AT > N . 

Since {7 is an open set in Ok and Ok is compact, t/ c fl Ok is a compact set in 
Ok- The relation 

f](0 K n O^o,cogj) c Ok n frgLr^coffL] c 17 

iV 
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holds. Taking their complement sets, we obtain 

\J(Ok n (frglrVgL, oo)) Da c no K 

iV 

Since K D (o-^]v) _1 ( C(7 2!L' 00 ) is °P en and ° K n (^}3v)~H co i2o> 00 ) C n 
(°i^]v+i) _1 ( ccr i^L; the compactness of ?7 C fl Ok guarantees the existence of an 
integer N such that 

Ok n 00) D f/ c n Ok 

for N > N . Thus, since <J^ N < crj^L, we obtain 

Ok n (4 n Jv) _1 [o, c^|Jv] c K n (ogkrMo, c^gj c C/ 

for AT > N . 

4.4 Proof of Theorem 12.61 

In order to show (I20|) . it is sufficient to prove that (c^L) c o a /?]v] contains the 
subspace V n £ 2 (Z+ ) . 

Since maxj e o ifn y n ^«+\ lla 7 !^ 2 is finite, we can choose c\ such that 

(^gL)~ 1 [o^iagL] = ^n£ 2 (z + ). 

Next, we fix an integer N and choose Co such that 

CqCt KNq > c x o K . 



For any N > N 



Thus, 



r n m > „ _(«) > „ _(«) 



Therefore, the set (er^J 1 [0, Cocr^jy] contains the subspace V fl £ 2 (Z + ). 
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Figure 2: Numerical result of functions for the ODE (I36j) 

5 Numerical examples 

Though the abstract structure of the algorithm is roughly explained in Subsection I2.2I 
of this paper, the detailed explanation of the practical algorithm requires many pages, 
which is reported in our paper [7], and so we omit it here. Here, we provide some 
numerical results only. In order that we can observe the accuracy of the algorithm, 
we chose example ODEs whose exact solutions are known analytically and can be 
written with special functions. Though the proposed method can be used widely 
even for the ODE's which can not be solved analytically, by intent we here give some 
examples for ODE's whose true exact solutions can be obtained analytically, in order 
to clarify how accurate solutions the proposed method gives. 
In the following, we use the bilinear form 

oo 
n=0 



32 



where 



w n : = 



1 (n<K) 
e r( M „-MK) (K < n < J) 

R := e'bJ-VK) ( n > ]V) 



with fi n : = 



fc + 1 A; + 1 



under the choice K = 2 



3(N-k ) 



J = 2[ 



15(Af-fc ) 
32 



J +fc , J = 2[ 



7(N-k ) 
16 



J + £; or if = 2 [ 



7(N-k ) 
16 



J +^0, 



J + ko, and r 



10 s 



The weight number series {w n }^L used in 



this bilinear form may seem to be somewhat complicated, but it is suitable for the 
symmetry property due to ifj ko ,h = ipko-n-k-i in fl2ZJ). 
The first example is the third-order ODE 

/'" - xf" - (81a; 2 - 54a; - 18u)f + (81s 3 - 54a; 2 - (18u + 162) x + 54)/ = .(36) 



is 



If v G Z + , the space of solutions in C 3 (R) fl L 

{C(exp ~ (3 7 1)2 )^(3x-l) | C G C}, where is a Hermite polynomial, because the 
differential operator on the left hand side of this ODE can be decomposed as 
^ x ) ' (l 5 (J^) 2 — (3a; — 1) 2 + {2u+ 1)) and it can be shown that there is no solution 
/ inL 2 fco) (R) such that ±f"+ (-(3a;-l) 2 +(2z/+l))/ belongs to ker(£-a;)\{0}. The 

results with v = 3, k = 4, N+l = 24, 36 , K = 2 + k and J = 2 [ m T f nl \ +h 

are shown in Figure [2], under the normalization (/, ^(ipk ,o + ^ko,-ko-i))n = 1- The 
errors of the result with N + 1 = 36 only are hardly noticeable in this figure. 



h 

Table 4: Numerical results for the ratio — under v = and k = 3 for the ODE (13 7p 

Jo 



ratio ^ 
io_ 



decimal expression of ratio 
iQ_ 



N+l 



50 



147826 
391819 



0.3772813467442875409308890074243464456802758 



100 



208588565 
552872013 



0.3772818303248061138518870912715200145245912 



150 



1969523740562 
5220298414229 



0.3772818303248061138245150519347658988268210 



200 



531796829098893 
1409547946268876 



0.3772818303248061138245150770765762118573286 



250 



651719569462020954 
1727407781341996633 



0.3772818303248061138245150770767548665927969 



300 



150649258697699321707 



399301653535776433703 



0.3772818303248061138245150770767548664028748 



true 



3+2v / 2CT(Erfc(^ 5 )-l) 



0.3772818303248061138245150770767548664028706 
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Number of significant digits Proportion p defined in (JM} 



10* 



10 1 



1 ♦ ♦♦♦♦♦ ♦ ♦ ♦ ♦ ♦ ♦♦ ♦ 



0.5- 



10* 



10 d 



10 4 
N+1 



10* 



10 J 



10 4 
N+1 



Figure 3: Number of significant digits of Figure 4: Proportion p denned in ([38 
the ratio for ODE (ETJ) for ODE (E 



Another example is Weber's differential equation (which is equivalent to the 
Schrodinger equation for a harmonic oscillator [3]) 

f'-a?f+(2v + l)f = 0. (37) 

As is well known, for v 6 Z + , the space of solutions in C 2 (M) fl L 2 (R) is 
{C(exp =§-)H v (x) | C e C}, which is a subspace of L 2 fco ^(IR) for any A; G Z + . For 
this example, convergence is very rapid, and we will report its accuracy by showing 
within how many digits the ratio between two coefficients /„ and /„/ in the expansion 
f( x ) = fn&n{x) coincides with the true ratio. For example, In Table IU we show 



h 

the results of the ratio — - for the case with v = 0, k = 3 K = 2[ 



J = 2[ 



15(iV-fco) 
32 



/o 



7(N-k ) 
16 



J + k and 



J + k Q , where the true ratio is obtained analytically (not numerically) 



by means of the computer algebra software package "Methematica" . Similar accuracy 
is observed for other ratios between the coefficients with small n and n' . With 
N + 1 = 7000, we obtained a result where it coincided with the true value up to 340 
digits. In Figure |3l we plot how the number of significant digits of this ratio depends 
on N. Moreover, we found that the rational ratios obtained in this case have almost 
a 'full precision', because the proportion 



P - 



(number of significant digits of the ratio) 



(number of digits of numerator) + (number of digits of denominator) 



(3* 



almost equals 1 for + 1 > 100 as is shown in Figure HI (In this case, the ratio j£ 
has no imaginary part due to a symmetry.) 



Number of significant digits 




N+1 



Proportion p defined in ( 138|) 



0.5 



♦ ♦ 



10^ 



10 J 



10 4 



N+1 
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Figure 5: Number of significant digits of the 

ratio s f and Number of significant digits of Fi S ure 6: Proportion p defined in p§, 
the ratio for ODE for ° DE ® 



Moreover, under the scale change x — > 30x, the accuracy is improved very much. 
Under this scale change, ODE (j3"Tj) is modified to 

^/»-(30)V/+(2i/+l)/ = 0. (39) 

The results for this ODE with v = 0, k = 6, K = 2[ 7{N ~ ko) \ + k and J = 
2|_ 15 ^~ fc °^ j + k are given in FigO In this case, with N > 100, the number of the 
significant digits between two coefficients is infinite (i.e. perfectly exact ratio is ob- 
tained) when the true ratio is rational which occurs for the ratio among f , fi, . . . , / 5 , 
and it is very large even when the true ratio is irrational. For example, for the ratio 
(which is irrational), the ratio obtained numerically by the proposed method co- 
incides within 8783 digits to the true ratio when N + 1 = 30000. Moreover, for the 
ratio between the values of the solution function at two points , there the nu- 

merical result by the proposed method coincides within 2599 digits to the true ratio. 
(The number of significant digits seems to be proportional to iV 0,79 empirically in 
the case when N is sufficiently large.) As for the ratio defined in (138!) . the numerical 
results by the proposed method give almost a 'full precision' al so in this case, which 
is shown in Fig EJ There results show evidently how accurate the proposed method 
is. 

Next, we give examples in Fuchsian class where Pm{%) has zero points. Here, 
we show how numerical results converge to true solutions in such a case where the 
true solutions are written with the associate Legendre function. In this case, we are 
successful to extract only the true solutions defined only in (—1,1) by the proposed 
method, where the obtained solutions is almost zero outside these intervals. 




An example of such a case is for the associate Legendre differential equation 

(1 _ x *)f> _ 2xf + (u(u + 1) - f = 0. (40) 

By means of the discussion in Subsection 13.5} we can treat this ODE by the proposed 
algorithm as the Fuchsian-type ODE 

(1 _ x y f » _ 2x (l - x 2 )f + [y{v + 1)(1 - x 2 ) -v 2 )f = (41) 

whose coefficient functions are polynomials. As is well known, there are three in- 
tervals (— oo, — 1), (—1,1) and (l,oo) within which smooth solutions are defined, 
because the coefficient function p2{x) of the highest order term has two zero points 
x — ±1. However, none of the solutions defined in the intervals (— oo, — 1) and 
(1, oo) is square-integrable, and hence the space of solutions in L 2 k JM.) (c L 2 (M)) 
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is the one- dimensional space {C ■ l[_i 5 i](x) • (1 — x 2 )^L^(x) \ C G C} indi- 
cator function, (1 — x 2 )^ L^(x): associate Legendre function). The results for this 
ODE with // = 3, v = 4, k = 6, K = 2[^^J + fc , J = 2L ^^ J + k 
and N + 1 = 50, 100 are given in Fig [71 Note that there the solutions are normal- 
ized by (/, o-(/0fc O) o + ^ko,-ko-i))n = 1- Surprisingly, almost only the component in 
{C ■ • (1 — x 2 )%L^(x) | C G C} is 'automatically' extracted, and the numer- 

ical solutions are almost zero outside the interval (—1,1), in spite of the existence 
of singularities at x = ±1. However, the convergence to the true solution is not so 
rapid as the cases where Pm( x ) has no zero points, though it converges to the true 
solution anyway. 

For the ODEs whose exact solutions can be written by the (associate) Laguerre 
functions x^e~^L^{x) within the interval (0, oo), we have already had similar results 
to this, where the obtained numerical solutions are almost zero for x < 0. 

6 Discussion 

6.1 Some properties of the basis functions used in this study 

The basis systems {e n }^ =0 and {e^}^L are closely related to Fourier series, by the 
change of variable 8 = 2 arctanx, as is shown in subsection 2.4 of the paper [8]. (The 
same change of variable has been used for a description of analytic unit quadrature 
signals with nonlinear phase [15] |16j.) 

The function ip k>0 is identical to the Cauchy wavelet [T7J used for continuous 
wavelet transformation [18J . Moreover, when k is even, ipk,n is closely related to the 
number state associated with su(l, 1) in a representation of su(l, 1) which can be 
formulated by adding a third generator to the two generators of the ax + b group |T9] . 

6.2 Extension to inhomogeneous differential equations 

The algorithm proposed in this paper is easily extended to linear inhomogeneous 
ordinary differential equations with inhomogeneous terms in T-fi . This extension only 
requires substitution of the right hand side of the simultaneous linear equations 
12 n bmfn = (m G Z + ) by the if-inner-products between the inhomogeneous term 
and the basis function e^. 

6.3 Modification of the method for the eigenvalue- eigenvector 
problem 

We have already proved that the proposed method can be applied for eigenfunction 
problems of self-adjoint operators with given eigenvalues, under some conditions, 
which will be reported in another paper [20] . In order to apply the proposed method 
to the eigenvalue-eigenvector problem for a linear operator, we must have a method to 
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obtain the eigenvalues, because the eigenvalue is regarded as a fixed parameter of the 
characteristic equation in the proposed method. In the case of discrete eigenvalues, 
if an eigenvalue is not exact, the function satisfying the characteristic equation does 
not belong to "H, and hence its corresponding vector is not square-summable. 

However, when we truncate the algorithm within a finite number of dimensions, 
the square-summability is not distinguishable. The number sequence obtained by our 
method for an approximate eigenvalue decays within a finite number of dimensions 
as rapidly as the number sequence corresponding to the true eigenvector. As the 
approximation of the eigenvalue is better, it decays for more dimensions. From this 
fact, we can propose a method to find the eigenvalue by observing the location of 
the bottom of the valley of the ratio Here we give an example of such valleys 

in Figure |HJ In this example, we are successful to separate two eigenvalues which 
are very contiguous by the 'tunnel effect', for a Schrodinger equation with quantum- 
double-well-type potential function. 




1738.8 



Figure 8: Example of 'valleys' of the ra- 
tio between the two norms 
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Figure 9: Possibility of a more precise in- 
terpolation by means of an index almost 
linear to the deviation 



Moreover, we have already invented another faster and more effective method 
for finding eigenvalues in a very high accuracy, based on a more analytical idea. 
This idea utilizes linear interpolations by means of some indices almost linear to the 
deviation of the eigenvalue (see Figj9]) which are calculated directly from numerical 
results. 



6.4 Possibility of the extension to partial differential equa- 
tions 

A similar idea to the proposed method can be applied to linear partial differential 
equations. However, the number of linearly independent solutions of simultaneous 
linear equations is not fixed but increasing as N increases for linear partial differential 
equations, while it is fixed at p = j + £ — 1 for linear ordinary differential equations. 



38 

Therefore, we have to estimate how much memory and how many calculations would 
be required. 

6.5 Possibility of the extension to weakly non-linear differ- 
ential equations 

This algorithm has the possibility of extension to nonlinear differential equations 
because of the following properties: ^From the definition of ip ko ,n, the relation 
V'feo,™i( x ) 'V^o,^ ( x ) = ^2*0+1, #ii+fi 2 ( x ) holds. The combination of this fact and Lemma 
13.21 results in the fact that the product i^k ,ni( x ) ' V'fco^O 3 ') can be expressed as 
a linear combination of i^ko,ni+n2 

(x),ip k0t n l+ n 2+ i(x), . . . , ip k0j n l+ii2+ko+1 (x). Similarly, 
the product of more than three basis functions can be written as a linear combination 
of finite numbers of the same basis functions. If the nonlinearity is weak, we can 
apply the proposed method to the successive approximation method for nonlinear 
differential equations, because of this property. However, for the nonlinear case, it is 
more difficult to find a proof of convergence and an upper bound for errors, than it 
is for the linear case. 

7 Conclusions 

We have proposed an integer-type algorithm which can determine accurately a basis 

M 

system for the space of solutions in T-L of the M-th order ODE Pm{x) (j^) m ) f{ x ) — 

m=0 

with polynomials or rational functions for the coefficient functions p m (m = 
0, . . . , M) under certain conditions. The basic structure of this algorithm has been 
shown in a more general framework and several conditions have been stated for the 
validity of this structure. Next, we have provided choices for the spaces and their 
basis systems satisfying these conditions, with detailed checks of these conditions. 
Thus, the validity of the proposed method has been proved. 

Moreover, we have shown convergence of the results of this method to true solu- 
tions of the differential equations, under the conditions required for the structure of 
the algorithm. Numerical results have indicated that this method has high accuracy. 
We have provided examples to show how the results converge to true solutions as 
the dimension of the subspace increases. 

This method will be extended or generalized for inhomogeneous equations, partial 
equations and weakly nonlinear equations in the near future, as has been mentioned 
in Section EJ Analyses of the accuracy and the amount of calculations required are 
also future problems. Moreover, it is our intent to apply this method, with some 
modifications, to the scattering problem in quantum mechanics. 
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A Proof of Lemma 3.1 



Proof of Lemma 13.11 ^From the last property of fl2T|) . { y ^ ipk,n | G Z} is or- 

thonormal. Therefore, we have only to prove the completeness in L? fc v(R). Let T be 
the Fourier transformation, where the Fourier transform of a function / is denoted 
1 f°° 

by \Tf)(y) '■— / f(x)e~ tyx dx. Some calculations by residue calculus result 



in 



iV2n e- y L, h (2y) (y > 0) 
V n>0, (.FVo,fi)(y) = { (42) 

(y<0) 

where L n (x) denotes the Laguerre polynomial of degree n. On the other hand, since 
V'o.ri^) = ^o, -n-i( x ) from (1271) . a property of the Fourier transform leads us to 

f -iV2^ e y Ln{-2y) (y < 0) 
V n>0, (^ ,-fl-i)(2/) = < (43) 

I (y > 0). 

Here, let 

-i 

£ (0) := ( E &^0,»(s) UnGC, fe}G£ 2 (Z\Z+)}, 
ri=— oo 
oo 

£ 5) : = { X;&V>o,ft(aO Un G C, {&} G £ 2 (Z+) } . 

ri=0 

Then, from the well-known fact that the set { e~2 L n (t), t > | n e Z + } is complete 
in L 2 ( M+ ) , we can show that { J 7 / | / G £+ 0) } = L 2 ( M+ ) . Similarly, from (gSJ and 

this fact, {J 7 / | / G £(o)} = £ 2 (K~). Since the null functions in L 2 (M) which are 
nonzero only at y — in the frequency domain belong to the kernel of the inverse 
Fourier transformation, from the Planchrel theorem, 

oo 

L 2 (R) = C © £+ = { MoM | G C, fe} G i\7L) } , (44) 
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and hence {y^^o,™ | n G Z} is complete in L 2 (IR) = L 2 ^(R). Then, since ^pk,h(x) = 
f 0, , X l , from (|23D and g4j), 

(x + i) K 

1 oo 

L (*)W = { T^ji E &V*,*0»0 I in e C, {&} e £ 2 (Z) } 

oo 

= { ^ Mk,n{*) I & € C, {&} G £ 2 (Z) } , 



and hence {y- ^fc,™ | n G Z} is complete in L 



2 I 



B Proof of Theorem 13.2 



For the proof of Theorem 13. 2\ here we start with the following lemma which is based 
on the translation of Lemma [3.21 by the 'matching' used in (l3"Tj) : 

Lemma B.l Let ko,j,m G Z + , k G Z and £1 := 2m + fco — k. Under the choices 
(1241) . ff25l) . and fl3~Tl) . /or k < k + m — j, the function x° (-^) m e n (x) can be expressed 
as a linear combination of e®, (n' = 0, 1, ...n + i\) at most for n < l\, and it can 
be expressed as a linear combination of e®> in' — n — t\, n — t\ + 2, n — i\ + 4, n — 
t\ + 6, ....,n + £i) for n > i\. In these linear combinations, all the coefficients are 
polynomials of iik ,n and ko with degree not greater than m. In particular, in the 
linear combination for n > l\, with n k0tn defined in (13~T1) . the coefficient of the first 

ko-K-j+m / -i \ j 



/ -\ Ko-K-j-i-m /jXJ m 

term with e^_ tl is ( — J ( - J (— l) m Y\ (hk , n + k + t) when n + k is even, 

and is ( — - J ( - ) II (^*o>« — * + 1) w/ien n + A; «s odd. 



The proof is derived directly from Lemma [3.21 together with (|3ip . 
Proof of Theorem I3.2t From the definition of Sq, the inequality k$ < k — s 
implies that v m G {0, 1, 2, M}, /cq < /c + m — degp m . Hence for every term in 

M deg p m 

the expansion P(x, J-^) = p m j x^ (-^) m , k$ < k + m — j holds. Therefore, 

m=0 j=0 

we can apply Lemma [B.ll term- wise in this expansion, where (x°(4-) m e r , e^)^o = 
for 

|r — n | > 2m + k — k^ and of course for |r — n| > 2M + fc — • H ence 
6; = (Be r , et) H o = for \r - n\ > 2M + k - k$ I.e. (a) holds. 

Next, we will show (b). Since n fco>r := [-^\ + (-l) r+ko+1 [^\, 

\rik ,r\ < — — ' ^ence, f° r fi xe( i ^0) f° r an y polynomial B(x), there exists 
a polynomial A(x) of the same degree as B(x) such that \B(hk , r )\ < Aw(?") for 
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r G Z + . Since Lemma [B.ll implies that there exists a polynomial Bi n ){x) of degree 
not greater than M such that h r n = B^(r) for every n G Z + , this fact results in the 
existence of a polynomial A(x) of degree not greater than M such that < A{r) 
for r G Z + , i.e. (b) holds. 

Moreover, for m < M — 1, ( xH4-) m e r , ,„„„ , , . ) =0 because 
\r - (r - (2M + A; - fc£))| > 2m + A; - fc£. Hence 

other hand, with n fc0jr := L-^J + (-l) r+feo+1 [^J, Lemma EH implies that 
T UA-\ M e P \ 

X W e ^ e r-(2M+k ^)/ n 



e r , e v . , , a, ) • On the 



k -k$+M 



M 



•1) M II ^ fe o,r + + t) (if fc + r : even) 



t=i 



k -k$+M Af 



t + i) 



(if fc + r : odd) 



l\3 



'±i V 

(=r) j ( — I . These facts and the 



holds for r > 2M + k — , because 

relation ^mj = Pm(±j) for r > 2M + fc — fco imply that 

i=o 



Be r , e 



'" "r-(2Af+fc -fc?)/ W 

fco-fe^+M 



-1) M JJ (n fco>r + fc + t) (if fc + r:even) 



t=i 



n 

t=l 



t + 1) 



(if k + r : odd). 



^From the definition of rifc 0ir , at least with r > ko + 2M, nfc 0ir + &o + ^ < ~2 for 
t = 0,1, 2, M when fc + r is even and rik ,r — t + 1 > 1 for t — 0, 1, 2, M 
when fc + r is °dd (where r = fc + 2M is impossible). Since pj^(±i) 7^ from the 

condition, we have the conclusion (Be T , , ,, , >6 , ) 7^ at least for 

r>k + 2M + max(-k$, 0) i.e. (c) holds. ■ 
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